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ABSTRACT 

The fuel optimal control problem arising in coplanar orbital transfer employing 
aeroassisted technology is addressed. The mission involves the transfer from high 
energy orbit (HEO) to low energy orbit (LEO) without plane change. The basic 
approach here is to employ a combination of propulsive maneuvers in space and 
aerodynamic maneuvers in the atmosphere. 

The basic sequence of events for the coplanar aeroassisted HEO to LEO orbit 
transfer consists of three phases. In the first phase, the transfer begins with a deorbit 
impulse at HEO which injects the vehicle into a elliptic transfer orbit with perigee inside 
the atmosphere. In the second phase, the vehicle is optimally controlled by lift and 
drag modulation to satisfy heating constraints and to exit the atmosphere with the 
desired flight path angle and velocity so that the apogee of the exit orbit is the altitude 
of the desired LEO. Finally, the second impulse is required to circularize the orbit at 
LEO. The performance index is maximum final mass. 

Simulation results show that the coplanar aerocapture is quite different from the 
case where orbital plane changes are made inside the atmosphere. In the latter case, 
the vehicle has to penetrate deeper into the atmosphere to perform the desired orbital 
plane change. For the coplanar case, the vehicle needs only to penetrate the 
atmosphere deep enough to reduce the exit velocity so the vehicle can be captured at 
the desired LEO. The peak heating rates are lower and the entry corridor is wider. 
From the thermal protection point of view, the coplanar transfer may be desirable. 
Parametric studies also show the maximum peak heating rates and the entry corridor 
width are functions of maximum lift coefficient. 

The problem is solved using a direct optimization technique which uses 
piecewise polynomial representation for the states and controls and collocation to 
represent the differential equations. This converts the optimal control problem into a 
nonlinear programming problem which is solved numerically by using a modified 
version of NPSOL. Solutions were obtained for the described problem for cases with 
and without heating constraints. The method appears to be more robust than other 
optimization methods. In addition, the method can handle complex dynamical 
constraints. 


* Staff Manager, and ** Senior Specialist, Advance Flight System, Advanced 
Technology. 
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NOMENCLATURE 


A : Sps/2 

Cp : drag coefficient 

Cdo - zero-lift drag coefficient 
Cl : lift coefficient 

Clr : lift coefficient for maximum lift-to-drag ratio 

D : drag force 

g : gravitational acceleration 

g s : gravitational acceleration at surface level 

H : altitude 

J : performance index 

K : induced drag factor 

L : lift force 

m : vehicle mass 

R : distance from Earth center to vehicle center of gravity 

Ra : radius of the atmospheric boundary 

Rc : radius of the low Earth orbit (LEO) 

Rd : radius of the high Earth orbit (HEO) 

Re : radius of Earth 

S : aerodynamic reference area 

t : time 

V : velocity 

T : thrust 

P : inverse atmospheric scale height 

Y : flight path angle 

¥ : heading angle 

a : bank angle 

0 : down range angle or longitude 

<j) : cross range angle or latitude 

P : gravitational constant of Earth 

p : density 

Ai : orbital plane changes 

AV : characteristic velocity 

subscripts 

c subscript for circularization or reorbit 

d : subscript for deorbit 

s : subscript for surface level 


1. INTRODUCTION 

In order to have a viable and affordable space program, advanced technology must be 
exploited and new design concepts must be developed to reduce the size and cost of 
transportation elements for supporting new mission requirements. One of the new 


462 


concepts that has evolved in recent years to advance the cost effectiveness of space 
ransportation systems is the aerodynamically assisted orbit transfer. Such an orbital 
ransfer vehicle is designed with an aerodynamic configuration which can utilize the 
planetary atmosphere for the purpose of energy management. Numerous studies 
lave demonstrated that the use of the aerobraking can significantly reduce the 
propulsive velocity requirements for certain class of orbit transfers. Excellent review 
papers were given by Warberg (Reference 1) and Mease (Reference 10). 

n our earlier studies, the fuel optimal control problem arising in a typical nonplanar 
prbit transfer from HEO to LEO as discussed in most recent publications was 
addressed. As discussed in References 2 and 15, the aeroassisted orbit transfer 
vehicle (AOTV) maneuver involves three phases with three propulsive burns or 
mpulses as sketched in Fig.1 . The orbital plane change was assumed to perform 
entirely inside the atmosphere with aeroassistance. Unlike References 2 and 15 , the 
nore general formulation given in Reference 1 7 does not restrict the orbital plane 
change to be performed entirely inside the atmosphere. In the first phase , the orbital 
ransfer begins with a deorbit impulse at HEO which injects the vehicle into an elliptic 
ransfer orbit with a plane change at HEO and with the perigee inside the atmosphere, 
n the second phase, the vehicle is inside the atmosphere and is optimally controlled 
py the lift and bank angle modulations to perform another orbital plane change and to 
satisfy the heating rate and other physical constraints. Because of the the energy loss 
luring the atmospheric maneuvers, an impulse is required to initiate the third phase to 
poost the vehicle back to the final orbital altitude. Finally, the third impulse is applied 
o circularize the orbit at LEO. Additional plane changes are allowed at the 
itmospheric exit and the final orbit circulation, in summary, there are three propulsive 
plane changes associated with three propulsive burns outside the atmosphere and an 
leroassisted orbital plane change inside the atmosphere. In Reference 17, simulation 
esults for the general formulation were obtained under the assumption that all 

rajectories enter the atmosphere at the same <j> e . Ve. and 0 e . In addition, simulation 
esults were compared with those obtained in Reference 2 and 15, where orbital plane 
vhanges are performed entirely inside the atmosphere. These studies provided 
lecessary data base and essential information concerning the effective use of 
peroassisted orbital plane changes. 

n this paper, the fuel optimal control problem arising in a typical coplanar 
leroassisted orbit transfer is addressed. The mission involves the transfer from high 
mergy orbit (HEO) to low enery orbit (LEO) without plane change, The basic 
approach here is to employ a combination of propulsive maneuvers in space and 
aerodynamic maneuvers inside the atmosphere. The aeroassisted orbital transfer 
aroblem is formulated under the assumption that no orbital plane change is needed. 
Similar to Reference 15 and 17, the basic sequence of events consists of three 
aliases but only two impulses are needed. In the first phase, the transfer begins with a 
jeorbit impulse at GEO which injects the vehicle into an elliptic transfer orbit with 
aerigee inside the atmosphere. In the second phase, the vehicle is optimally 
aontrolled by lift and drag modulation to satisfy heating and other physical constraints 
and to exit the atmosphere with the desired flight path angle and velocity so that the 
apogee of the exit orbit is the altitude of the desired LEO. Finally, in the third phase, 
he second impulse is required to circular the orbit at LEO. The optimal control 
solutions were all obtained by using the Hermite polynomial and collocation technique 
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to convert the optimal control problem into a corresponding nonlinear programming ( 
NP ) problem which is solved numerically using the optimization code, NZSOL ( cf. 
Reference 12 ) provided by Gill, which is an improved version of NPSOL ( cf. 
Reference 6 ), developed at Stanford. This solution method is different from the 
indirect method such as those discussed in Reference 2,4,7 and 8. Simulation results 
were then compared with those obtained earlier for different orbital inclination 
changes in Reference 15 and 17. The details are presented and discussed here. In 
this paper, simulation results were actually obtained for returns from geosynchronous 
orbit (GEO) to space station orbit (SSO). It is important that in the future these 
simulations be extended to include all other realistic flight constraints and to establish 
baseline optimum trajectory characteristics for GEO to Space station or shuttle lunar 
and Mars missions. 

2. DIRECT TRAJECTORY OPTIMIZATION WITH COLLOCATION AND 
HERMITE POLYNOMIALS 

In the direct collocation with nonlinear programming approach, the trajectory is 
approximated by piecewise polynomials, which represent the state and control 
variables at a number of discrete time points, i.e., nodes. For a given state variable, 
the state trajectory over a given "segment" between two nodes is taken to be the 
unique Hermite cubic which goes through the end points of the segments with the 
appropriate derivatives that are dictated by the differential equations of motion at the 
endpoints. This is the "Hermite cubic" since it is determined by the states and their 
derivatives. A collocation is taken at the center of the segment where the derivative 
given by the Hermite cubic is compared to the derivative obtained from the evaluation 
of the equations of motion. The difference is termed the "defect" and is a measure of 
how well the equations of motion are satisfied over the segments. If all the defects are 
zero, then the differential equations are satisfied at the center collocation points as 
well as at the endpoints. Figure 2 shows the typical defects between node 1 and node 

Let the system of equations of motion be given as 

X' = f(X,U,D) (2-1 a) 

where X is the state vector, U is the control vector, D is the design parameter vector 
and (’) denotes the differentiation with respect to the time. Let the time over a given 
segment be T. For the problems mentioned above, one can show that 
X = ( x, y, z, x’, y’, z’, m ) 

u = ( c L,o) (2-lb) 

D = (Ai-j, Ai' 2 , A^) 

where design parameters are defined here as unknown constants ( i.e., three 
propulsive plane changes ) to be determined by the optimizaton processes. Then the 
Hermite interpolated x-component of the state vector X at the center point is 

x c = (1/2) ( X1 + x r ) + (T/8) [f(X t ,Ui) - f(X r ,U r )] (2-2) 

where xi and x r are respectively the x-component of the state vector X at the left and 
the right nodes. The derivative of the interpolating Hermite cubic at the center point is 
x c ' = -3/(2T) ( Xl - x r ) - (1/4) [f(X 1( Ui)+ f(X r ,U r )] (2-3) 

The defect vector is then calculated as 

d = f(X c Uc) - x c ' (2-4) 

If *1, ui> Xr> and Ur are chosen such that the elements of the defect vector, d, are 
sufficiently small, the "Hermite polynomials" become an accurate approximation to the 
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solution of the differential equations of motion (by implicit integration). With the above 
approach, the differential equations are converted into nonlinear algebraic constraint 
equations and the optimal control problem can then be solved using the nonlinear 
programming techques. 

3. BASIC EQUATIONS FOR OPTIMAL AEROASSISTED ORBITAL 
TRANSFER 

The aeroassited orbital transfer can be analyzed in three phases , i.e., deorbit, 
aeroassist (or atmospheric flight), boost and reorbit (or circularization). In each of the 
phases, a particular set of equations of motion apply. 

3.1 Deorbit 

Initially, the spacecraft is moving with a circular velocity V d = ^l7r^ in a circular orbit 
of radius Rd, well outside the Earth’s atmosphere. Deorbit is accomplished at point D 
by means of an impulse AVd, to transfer the vehicle from a circular orbit to an elliptic 
orbit and with perigee low enough for the trajectory to intersect the dense part of the 
atmosphere. Since the elliptic velocity at D is less than the circular velocity at D, the 

impulse AVd is executed so as to oppose the circular Velocity Vd ■ In other words, at 
point D, the velocity required to put the vehicle into elliptic orbit is less than the velocity 
required to maintain it in circular orbit. The deorbit impulse AVd causes the vehicle to 
enter the atmosphere at radius R a with a velocity V e and flight path angle y e . It is 
known that the optimal energy loss maneuver from the circular orbit is simply the 
Hohmann transfer and the impulse is parallel and opposite to the instantaneous 
velocity vector. 

After applying the deorbit impulse and before entering the atmosphere at R a , the 
deorbit trajectory is a coasting arc and known integrals of the equations of motion can 
be used to relate the state vectors at R a ,the entry into atmosphere to the state vectors 
right after the deorbit impulse at Rd. Using the principle of conservation of energy and 
angular momentum at the deorbit point D and the atmospheric entry point E, 

V e 2 /2-p/R a =V 1 2 /2-p/R d 

R a V e cos (-Y e ) = R d ^ 

where Vi is the magnetude of the velovity right after the deorbit impulse AV d 
the above equations we can solve for Vi and then compute AV d to get 

V, = ^2n(1 / R a - 1 / R d ) / [(R a / R a ) 2 / COS 2 Ye - 1] 

and 

AV d = V d -V 1 

It is easily seen that the minimum deorbit impulse AVdm obtained for y e = 0, 
corresponds to an ideal transfer with the space vehicle grazing the atmospheric 
boundary. To ensure proper atmospheric entry, the deorbit impulse AVd must be 
higher than the following minimum deorbit impulse AVdm 


we get 

(3-1) 
(3-2) 
and from 

(3-3a) 

(3-3b) 
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Vim = few / R a - 1 / R<J> / [(Rd / R a f - 1] 

^Vd m = V(j — Vi m (3_4b) 

Physically, the second term of the above equation corresponds to the apogee velocity 
of an ellpitic transfer orbit with perigee radius R a and apogee radius Rd This elliptic 
transfer orbit is tangent to the atmosphere boundary at perigee. It will be shown later 
that the nonlinear constraint equations ( 3-15 ) at the atmospheric entry point can also 
be derived from equations ( 3-1 and 2 ). 


3.2 Aeroassist 

During the atmospheric flight, the vehicle can be optimally controlled by the lift and 
bank angle modulations to achieve the necessary velocity reduction (due to the 
atmospheric drag) and the orbital plane change if needed. In the present formulation, 
only the aeroassisted atmospheric flight need be solved by using the collocation and 
nonlinear programming techniques discussed earlier in this paper. The solutions in 
the other phases are provided by the known integral relations of the equations of 
motion because these arcs are coasting arcs. 


Consider a vehicle with the point mass m, moving about a rotating spherical planet. 
The atmosphere surrounding the planet is assumed to be at rest, and the central 
gravitational field obeys the usual inverse square law. The equations of motion for the 
vehicle are given by (Figure 1), 
r = Vsiny 

^ _ Vcosycosy 
rcos<|) 

■ _ V cosy siny 


^ _ (t)T cose-D) 

m r 

. _ (r)T sine + L)coso 
y ~ mV 


-- l n ^ + w 2 r cos<J) (siny cos<}>- cosy siny sin<f>) 


ficosy Vcosy 
Vr 2 + r 


+ 2co cosij / cos<{> 


(3-5a) 


(3-5b) 


(3-5c) 


(3-5d) 


+ 

¥ = 


co r cos <|> 


(cosy cos<J) + siny siny sin<J>) 


(r|T sine + L)sina V cosy cosy tan<J> 
mV cosy r 


2co(tany siny cos<J> — sin<})) 


(3-5e) 


, co 2 r cos y sine}) cos<|) 
V cosy 
m = -f(r,V,ri) 


(3-5f) 

(3-5g) 


where for a given vehicle, the drag D and the lift L are 

D = ApV 2 C D 


(3-5h) 
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L = ^-pV‘ 
2m 


(3-Si) 

(3-53) 


and the drag and lift coefficients obey the drag-polar relation 

Cq = C D0 + KCl 
A lso, for an exponential atmosphere, one has 

P = Ps ex P (~ H P) and H = R - R e (3-5k) 

Simulation results obtained here were using the U.S. standard Atmosphere 1976. 


For aeroassisted orbital transfer problems considered here, one assumes that, inside 
the atmosphere, the vehicle is optimally controlled by the aerodynamic forces only. It 
is assumed that the thrust T is absent and the point mass is constant in this region. 
Furthermore, no earth rotation was assumed. The later is equivalent to consider the 
motion with respect to an earth fixed inertial coordinate system ( ECI ). The plane 

change or the orbit inclination, /', is related to the cross range <|> and the heading angle 

\jr as 

cosi = cos <|> cosy t e <t<tf (3_0) 

For coplanar orbital transfer problems considered here, the orbit inclination is 
assumed to be constant throughout the atmospheric flight. 


3.3 Boost and Reorbit 

During the atmospheric flight, the vehicle is optimally controlled by the lift and drag 
modulation to satisfy the heating constraints and to exit the atmosphere with the 
desired flight path angle and velocity so that the apogee of the exit orbit is the altitude 
of the desired LEO. Thus, no impulse is required at the exit from the atmosphere to 
boost the vehicle back to the final orbital altitude at LEO. The vehicle exits the 

atmosphere at point F, with a velocity Vf and the flight path angle Yf. The additional 

impulse AVt>, required at the exit point F for boosting the vehicle into an elliptic 
transfer orbit with apogee radius R c is assumed to be zero and the reorbit (or 
circularization) impulse AV C , required to insert the vehicle into a circular orbit are 
obtained by using the principle of conservation of energy and angular momentum at 
the exit point Fand the reorbit or circularization point C. Thus, we have 

V, 2 /2-n/R a = V|/2-n/R c (3 .7) 

Vf R a cosyi =R c V 3 ( 3 -8) 

where Vf is the velocity at the exit from the atmosphere and V3 is the velocity at the 
reorbit point C just before the circularization burn A V c . 


Solving for Vf and V 3 from the above equations (3-7) and (3-8) yields 


J 2p(1/R a -1/R c )/ 

1 -(R a / R c ) 2 cos 2 y f ] = 

^2n(1/R a -1/R c )/ 

(R c /R a ) 2 /cos 2 ff-1 


and AV C can be computed as follows 
AV b = 0 


(3-9) 

(3-10) 

(3-11) 
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av c = v c -v 3 


(3-12) 

It is interesting to note that V 3 is maximum for y? =0 and therefore the reorbit impulse 

A V c ,is minimum for yf =0. It will be shown later that boundary conditions and nonlinear 
constraint equations at the exit point F, can be derived in terms of the final orbit 
characteristics and the final state vectors at the exit as shown in (3-16,17,& 18). 

3.4 Performance Index 

It is known that the change in speed, AV, also called the characteristic velocity, is a 
convenient parameter to measure the fuel consumption. For minimum-fuel maneuver, 
the objective is then to minimize the total characteristic velocity. A convenient 
performance index is the sum of the characteristic velocities for deorbit, boost, and 
reorbit, as 

J = AV d + AVq 

= AV d (R d , Ah, V e , y e ) + AV C (R c , Ai 3 , y f , V f ) ^ 

Where, AVd and AV C are the deorbit, and reorbit characteristic velocities respectively, 
and are given by (3-3, and 12) respectively. Note that for a given final circular orbit, 
the impulse AV C are completely determined by the state variables Vf and yt at the exit 
of the atmospheric portion of the trajectory. The velocity V e and the flight path angle 
y e at the atmospheric entry point are dependent only on the magnitude of the deorbit 
impulse AVd. It follows that the optimal control problem needs to consider only the 
trajectory segment within the atmosphere subject to the nonlinear constraints and 
boundary conditions at the atmospheric entry and exit points. In addition, other path 
constraints such as the peak heating rate have to be satisfied. 

3.5 Boundary conditions and constraints 

The boundary conditions and constraints for the optimal control problem can be 
summarized as follows: 


At the entry into atmosphere, the following initial constraints must be satisfied. 


R = R a 


Ye - 0 


vi 


= 0, \j/ e = 0, e e = 0, 


cos 2 (y e ) 



(3-1 4a) 
(3- 14b) 


-w 


V^a 


_ 1 _ 

Rd 


= 0 


(3-15) 

The first initial constraint is required to ensure the vehicle enters the atmosphere. The 
second set of boundary conditions assumes that all trajectories enter the atmosphere 


at the same cp e > Ve and e e- In the present formulation, the initial velocity and the flight 
path angle are unknown and to be determined by the optimization processes subject 
to the constraint equation (3-15). 


• At the exit from atmosphere, the following constraints must be satisfied. 


R = R a ; y f > 0 


(3-16) 
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v, 2/ 


i-5i 


cos 2 Yf 


\ 

1 1 


M* 

) 

^Ra Rc ; 


= 0 


(3-17) 

cos if =cos<J>f cosy f = cosi e (3-18) 

Equation (3-16) is required to ensure the vehicle exit the atmophere. The second 
costraint (3-17) must be imposed to determine the correct Vf and yt if AVt>,\s assumed 
to be zero as in the coplanar case discussed here. The third constraint (3-18) is 
required to ensure the orbital transfer is coplanar. 


In addition, there are other path constraints ,i.e., constraints must be satisfied along the 
trajectory such as stagnation point heating rate constraints, altitude constraints, 
bounds on the control variables abd others 


4. STRUCTURE AND SOLUTION OF THE NONLINEAR PROGRAMMING 
PROBLEM 

The direct collocation and Hermite polynomial procedures described above convert 
optimal control problems into corresponding nonlinear programming problems. 
Ordinary differential equations are converted into corresponding nonlinear algebraic 
equations (or nonlinear “defects” constraint equations). These problems can then be 
solved using nonlinear programming codes. 

The variables for the nonlinear programming problem are the collected state vectors 
and control vectors at the nodes and the time duration of phases. These quantities are 
assembled into the NLP state vectors 

X T =[xW, xJ.ltf.M2 tk] ( 4. 1} 

where n is the number of nodes and k is the number of phases on the trajectory. The 
defects and other physical and mathematical constraints are collected into the NLP 
constraint vector C 

C T = [d|,d l, dj, w|, wj, wj, wT] 

where dj is the defect vector and w is a vector of additional problem constraints. 

The nonlinear programming code used here is the NZSOL (Reference 12). The 
NZSOL is an improved version of the NPSOL (Reference 6) , developed by the 
Stanford Optimization Laboratory and designed to minimize a smooth nonlinear 
function subject to a set of constraints which may include simple bounds on the 
variables, linear constraints, and smooth nonlinear constraints The problem is 
assumed to be stated in the following form: 

NP 

minimize F(x) 
xeR n 


x 

subject to ^<Mi_x>< u, 

l c (x)J (4-3) 

where the objective function F (z) is a nonlinear function, Al is an mu x n constant 
matrix of general linear constraints, and c(x) is an m/v - vector of nonlinear constraint 
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functions, the objective function F and the constraint functions are assumed to be 
smooth, i.e., at least twice-continuously differentiable. (The method of NPSOL will 
usually solve NP if there are only isolated discontinuities away from the solution). 

Note that upper and lower bounds are specified for all the variables and for all the 
constraints. This form allows full generality in specifying other types of constraints. In 
particular, the i-th constraint may be defined as an equality by setting = uj. If certain 
bounds are not present, the associated elements of £ or u can be set to special values 
that will be treated as - 00 or + 00 . 


Here we briefly summarize the main features of the method of NZSOL and NPSOL as 
discussed in Reference 6 because Reference 12 is not available to general public. At 
a solution of NP, some of the constraints will be active, i.e., satisfied exactly. An active 
simple bound constraint implies that the corresponding variable is fixed at its bound , 
and hence the variables are partitioned into fixed and free variables. Let C denote 
the m X n matrix of gradients of the active general linear and nonlinear constraints. 
The number of fixed variables will be denoted by hfx, with n fr (hfr = n - rrx) the 
number of free variables. The subscripts “FX” and “FR” on a vector or matrix will 
denote the vector or matrix composed of the components corresponding to fixed or 
free variables. The details are discussed in Reference 11. 


A point x is a first-order Kuhn-Tucker point for NP if the following conditions hold: 


(i) 

(ii) 



x is feasible; 


there exist vectors £ and X (the Lagrange multiplier vectors for the 
bound and general constraints) such that 

g = c T x+c, 


(4-4a) 


where g is the gradient of F evaluated at x, and Q= 0 if the j-th variable 
is free. 

The Lagrange multiplier corresponding to an inequality constraint 
active at its lower bound must be non-negative, and non-positive for 
an inequality constraint active at its upper bound. 


Let Z denote a matrix whose columns form a basis for the set of vectors orthogonal to 
the rows of Cfr\ i.e., CfrZ = 0. An equivalent statement of the condition in terms of Z 
is 

zT 9fr=° (4-4b) 

The vector Z r gFR is termed the projected gradient of F at x. Certain additional 
conditions must be satisfied in order for a first-order Kuhn-Tucker point to be a solution 
of NP. 


4.1 The Quadratic Programming Subproblem 

Similar to NPSOL, the basic structure of NZSOL involves major and minor iterations. 
The major iterations generate a sequence of iterates ( x ^ that converge to x\ a first- 
order Kuhn-Tucker point of NP. At a typical major iteration, the new iterate x is 
defined by 
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x - x + ap, (4-5a) 

where x is the current iterate, the non-negative scalar a is the step length, and p is the 
search direction. Also associated with each major iteration are estimates of the 
Lagrange multipliers and a prediction of the active set. 


The search direction p is the solution of a quadratic programming subproblem of the 
form 


minimize 

P 


T 1 Tii 

0 P+gP H P 


subject to i< 


P 

, a lP 

I a npJ 


SU, 


(4-5b) 

where g is the gradient of F at x, the matrix H is a positive-definite quasi-Newton 
approximation to the Hessian of the Lagrangian function and An is the Jacobian 
matrix of c evaluated at x. 


The estimated Lagrange multipliers at each major iteration are the Lagrange 
multipliers from the subproblem (and similarly for the predicted active set) and provide 
information about the the sensitivity of these NLP problems. 


Certain matrices associated with the QP subproblem are relevant in the major 
iterations. Let the subscripts "FX" and “FR" refer to the predicted fixed and free 
variables, and let C denote the m x n matrix of gradients of the general linear and 
nonlinear constraints in the predicted active set. First, we have available the TQ 
factorization (Reference 1 1 ) of Cfr : 

c fr Qfr = (0 T ). ( 4-6) 

where T is a nonsingular m x m reverse-triangular matrix (i.e., f,y = 0 if i +j< m ), and 
the non-singular pfr x pfr matrix Qfr is the product of orthogonal transformations. 
Second, we have the upper-triangular Cholesky factor R of the transformed and re- 
ordered Hessian matrix 


R t R = H q = Q t HQ, (4 _ 7) 

where H is the Hessian H with rows and columns permuted so that the free variables 
are first, and Q is the n x n matrix 

f Q FR ^ 


Q = 


'fx. 


(4-8) 


with Ifx the identity matrix of order npx- If the columns of Qfr are partitioned so that 

Qfr = ( z Y), (4-9) 


the n z (n z = p F r - m) columns of Z form a basis for the null space of Cfr ■ The matrix 
Z is used to compute the projected gradient Z T gFR at the current iterate. 


As discussed in Reference 6 and 1 1 , a theoretical characteristic of SQP methods is 
that the predicted active set from the QP subproblem is identical to the correct active 
set in a neighborhood of x*. In NPSOL, this feature is exploited by using the QP active 
set from the previous iteration as a prediction of the active set for the next QP 


subproblem, which leads in practice to optimality of the subproblems in only one 
iteration as the solution is approached. Separate treatment of bound and linear 
constraints in NPSOL also saves computation in factorizing Cfr and Hq. 


4.2 The merit function 

Detailed discussions of the merit function are given in Reference 14. In NZSOL and 
NPSOL, once the search direction p has been computed, the major iteration proceeds 

by determining a steplength a that produces a “sufficient decrease” in the augmented 
Lagrangian merit function 


L(x,X,s) = F(x)-£ Xi(Cj(x)-Si) + -I Pj (Ci(x)-Sj) 2 , 

1 1 (4-10) 

where x, X and s vary during the line search. The summation terms involve nnlv the 

QPQlinear constraint?. The vector X is an estimate of the Lagrange multipliers for the 


nonlinear constraints of NP. Ihe non-negative slack variable {s,} allow nonlinear 
inequality constraints to be treated without introducing discontinuities. The solution of 
the QP subproblem (4-5) provides a vector triple that serves as a direction of search for 
the three sets of variables. 


4.3 The quasi-Newton updated 

Before going into the detailed discussions, it is important to point out that both the 
NZSOL and NPSOL start by initializing the Hessian matrix H = Identity matrix. Thus at 
the beginning, the search direction is in the steepest decent direction. No initial 
curvature information is computed and the curvature information is accumulated 
through the BFGS quasi-Newton updates. The matrix H in (4-5) is a positive-definite 
quasi-Newton approximation to the Hessian of the Lagrangian function. At the end of 
each major iteration, a new Hessian approximation H is defined as a rank-two 
modification of H. In NPSOL the BFGS quasi-Newton update is used: 

H - H — — Hss T H + -^yy T , 

sHs * (4-11, 

where s = x-x (the change in x ). 

Rather than modifying H itself, the Cholesky factor of the transformed Hessian Hq (4- 
7) is updated, where Q is the matrix from (4-8) associated with the active set of the QP 
subproblem. The update (4-11) is equivalent to the following update to Hq : 

Hq = Hq - ~r u - HqSqSqH q + y Q yQ, 

s Q h Q s Q yo s Q ( 4 .12) 

where y Q = Q T y and s q =Q t s. This upda te mav be expressed as a mnk-nnp 

update tQ Q and is used to incorporate new curvature information obtained in the 
move from x to x . 


4.4 NZSOL, NPSOL 4.02, and NPSOL 2.1 

For those who are interested in applying these NLP codes, there are two publised 
versions of NPSOL. The NPSOL 4.02 was developed after the NPSOL 2.1 and 
therefore more reliable and efficient algorithm were incorporated according to Gill ( 
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Reference 12 ). However, in updating the Cholesky factor, the NPSOL 4.02 updates 
the whole or complete R while the NPSOL 2.1 updates only the part associated with 
the Z-space or null space of R. For the problem formulated here .usually several 
hundred varibles are involved and the NPSOL 2.1 converges in less computing time. 
The NZSOL (Reference 12) incorporates not only latest efficient and reliable algorithm 
but also updates only the part of R associated with the null space of R only. In addition 
to improve the algorithm of NPSOL, it also adopts the best parts of both NPSOL 2.1 
and 4.02. 

Finally, it may be interesting to point out that the matrices in the present formulation 
using collocation and Hermite polynomial are large and fairly sparse. For 
computational efficiency, it is important to incorporate NLP codes such as MINOS 
(Reference 13) to take advantage of the special characteristic of the collocation 
formulation discussed here. 

5. NUMERICAL RESULTS AND DATA 

The data used in the numerical experiments presented here (c.f. Reference 2 and 9) 
are summarized as follows: 

Cqo = 0.1 ; K = 1.111 ; m/S = 300 kg/m 2 (5-1) 

and the drag polar is 

C D = C D o + K * C L 2 ( 5-2 ) 

and other useful data are 

p a = 1.225 kg / m 3 ; ji = 3.986x1 0 14 m 3 / sec 2 

(3 = 1/6900 m" 1 ; R E =6378 km 

H a = 1 20 km (5-3) 

Using the above mentioned data, simulations were carried out 

For an AOTV returning from the geosynchronous orbit (GEO) to the space station orbit 
(SSO) , one has Rd = 42240 km and R c = 6934 km. Simulation results were obtained 
for the following parametric studies for different values of Clm . 

a) Case 1 ( Reference Case ). For this reference case , simulation results were 
obtained under the general formulation that no orbital plane changes are allowed at 
deorbit, boost, and reorbit impulses and inside the atmosphere. This reference case 
has the following entry and exit status 

Entry status: H e = 1 20 km; V e = 1 0.31 5 km/sec 

y e = -3.727 degrees; <J> e = 0; \jr e = 0 (5-4) 

Exit status: Hf = 1 20 km; Vf = 7.952 km/sec 

y f = 0.91 deg; <}>f = 0 deg 

= 0 deg; total flight time = 769.25 sec ( 5-5 ) 
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The characteristic velocities are 1489 meters per seconds and 132 meters per second 
at the deorbit and the reorbit respectively. The total characteristic velocity is 1621 
meters per second. The Ci_Max is assumed to be 3.0. 

b) Case 2 : Simulation results were also obtained for the case where the maximum lift 
coefficient is assumed to be 4.5. Similar to the reference case, no orbital plane 
change is allowed. The entry and exit status are summarized as follows:. 


Entry status: 

H e = 120 km; V e = 10.3118 km/sec 
Ye = -3.576 degrees; be = 0; \jr e = 0 

(5-6) 

Exit status: 

H f = 1 20 km ; V f = 7.951 8 km/sec 
Yf = 3.576 deg; bf = 0deg 

\|/ f = 0 deg; if = 0 degree , total flight time = 673.38 sec 

(5-7) 


The deorbit characteristic velocity is 1488.91 meters per second and the recirculation 
characteristic velocity is 131.61 meters per second. The total characteristic velocity is 
1620.517 meters per second. It is interesting to observe that all the characteristic 
velocities are almost the same as those obtained in case 1. 

c) Case 3 : Similar to Case 2 , the optimal control solution has a maximum lift 
coefficient of 2.3 and has the following entry and exit status. 

Entry status: H e = 120 km; V e = 10.31 15 km/sec 

y e = -3.814 degrees; <J> e = 0; \|/ e = 0 (5-8) 

Exit status: Hf = 1 20 km; Vf = 7.951 5 km/sec 

y f = 0.917 deg; <j>f = 0.deg; if =0.0 deg 

\j/ f = 0.0deg; total flight time = 857.34 sec (5-9) 

d) Case 4. Similar to Case 2, the optimal control solution has a maximum lift 
coefficient of 0.9 and has the following entry and exit status. 

Entry status: H e = 120 km; V e = 10.31 17 km/sec 

Ye = -4.154 degrees; <J> e = 0; \|/ e = 0 (5-10) 

Exit status: Hf=120km; Vf = 7.9509 km/sec 

Y f = 0.954 deg; <>f = O.deg; if = 0.0 deg 

\|/ f = 0.0deg; total flight time = 1450.67 sec (5-11) 

Again, all the characteristic velocities associated with the deorbit, and reorbit impulses 
for Case 3 and Case 4 are almost the same as Case 1 and Case 2. 

Time histories of altitude, velocity, flight path angles, ift coefficient, lift to drag ratio, 
dynamical pressure, atmospheric density and heating rate for all three cases ( Case 1 , 
2, & 3.) are shown in Figure 3-10. It is important to point out that for the vehicle 
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considered here, the maximum lift coefficient is 0.9. However, for the simulation 
results shown here, the maximum lift coeffient is set to be less than 4.5 for paramatric 
studies. 

An interesting observation from the simulation results as shown in Figs. 3 to 10 is that 
although the total characteric velocity is insensitive to the variation of the maximum lift 
coefficient , the optimal trajectory is very sensitive. The higher the maximum lift 
coefficient is the less the vehicle penetrates into the atmosphere and the less the time 
of flight is needed. The vehicle was flying at the maximum lift coefficient in contrast to 
previous simulation where the vehicle was flyiny at the maximum lift to drag ratio as 
shown in Reference 15 and 17. Simulation results shows that for the coplanar orbital 
transfer the vehicle has only to penetrate the atmosphere deep enough to reduce the 
exit velocity so that the vehicle will exit the atmosphere at the desired velocity and 
filght path angle so that the apogee of the exit orbit is the altitude of the desired final 
LEO. 

The heating rate Q r , along the atmospheric trajectory, is computed for the stagnation 
point of a sphere of radius of one meter, according to the following relation 
(Reference 2 and 6) 

Qr= K r p°- 5 V308 ( 5-12 ) 

where the p is the atospheric density in kg/km 3 , V is the velocity in km/sec and the K r is 
the proportionality constant equal to 0.000308. The time history of heating rates for the 
reference case ( Case 1 ) , Case 2, Case 3, and Case 4 were shown in Fig. 6. These 
simulation results presented provided enough information that the peak heating rate 
for coplanar case will be much less than those for cases aeroassisted orbital plane 
changes are made. As shown in References 15 and 17, one needs less thermal 
protection materials and more fuel consumption to fly the heat constrained trajectories 
and therefore by taking into account the weight of thermal protection materials one 
may find an optimal design to minimize the total vehicle weight. 

Another interesting observation from previous simulation results ( cf. Reference 15 and 
17 ) is that for given HEO and LEO, the deorbit impulse is almost the same for all the 
cases simulated here. The total characteristic velocity for a given optimal trajectory is 
almost completely determined by the boost and the recirculation. In fact, the boost 
velocity contributes the most to the variation of the total characteristic velocity. 
Physically, it is obvious as the vehicle makes a larger turn it also loses more energy 
and therefore needs more velocity to boost it back to the final orbital altitude. Although 
the total characteristic velocity is insentive to the magnitude of deorbit impulse, the 
optimal trajectory is very sensitive to AV(j . In the coplanar aeroassisted orbit transfer 
here, the boost impulse is not needed and the deorbit and recirculation impusles are 
almost the same for all the cases simulated here. Thus the total characteristic velocity 
is not sensitive to the variation of the maximum lift coefficient. However, the depth of 
penetration was shown as a function of the maximum lift coefficients. 

6. CONCLUDING REMARKS 

An excellent survey of the subject was given in Reference 1 . Walberg reviewed the 
problem of optimal aeroassisted orbital transfer with plane change. In a recent paper 
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by Naidu (c.f. Reference 2), fual optimal trajectories of aeroassisted orbital transfer with 
plane change were presented using the so-called multiple shooting method for the 
case without heating rate constraints and under the assumption that all the plane 
change was performed entirely in the atmosphere. A brief review of the progress 
made in this field was also given in Reference 2. In Reference 15 , a similar problem 
for cases with and without peak stagnation point heating rate constraints was solved 
using the coljocation and nonlinear programming technique. This method is 

especially suitable for parametrical studies because of its relative insensitivity to initial 
guesses. 


In Reference 17 , simulation results were obtained under a more general formulation 
that not all the orbital plane changes are made in the atmosphere. It must be noted 
that the AOTV transfer can be made more efficient propulsively if the plane change is 
performed partly in the atmosphere and partly in space and the propulsive plane 
change in space is subdivided into components associated with various impulsive 
points. All these plane changes were automatically determined by the optimization 
processes discussed in Reference 17 


The above studies provided necessary data bases and essential information 
concerning how to use and how to combine the propulsive and aeroassited orbital 
plane changes effectively. In this paper, another group of problems under the 
assumption no orbital plane change is allowed are investigated. In fact, the present 
investigation is closely related to the problem of returning from GEO to space station 
assuming all plane changes are made propulsively outside the atmosphere. As 
discussed, the charateristics of the flight is quite different from the cases where the 
orbital plane changes are made inside the atmosphere. In the latter case, the vehicle 
has to penetrate deeper into the atmosphere to perform the the desired orbital plane 
change. On the other hand, for the coplanar case, the vehicle needs only to penetrate 
the atmosphere deep enough to reduce the exit velocity so that the vehicle cn be 5 

capture at the desired LEO. 

It should be mentioned that the collocation and nonlinear programming technique 
discussed here was recently applied to another group of orbital transfer problem by 
Enright and Conway in Reference 3 and the relative insensitivity of this method to the 
initial guesses was also observed by them. Our basic simulation test bed is the OTIS 
fjj es j Reference 5 ) with an improved and updated nonlinear programming code 
(NZSOL ). All physical rnodals usGd wgtg docurnGntGd in RGfarGncG 5. Of cours© 
necessary modifications and corrections have to be incorporated to simulate the 
aerobraking problems discussed here. 

It may be worthwhile mentioning that the present problem was actually solved by 
guessing the initial state and control varibles at four selected points , i.e., the initial 
point, the final point and two other nodal points along the trajectory inside the 
atmosphere. The initial state and control variables at other nodes or grid points were 
simply obtained by linear interpolation. These initial guesses do not have to satisfy 
either the governing equations or the nonlinear constriants including the defects. Only i 
rough guesses are needed at these four points. Converged solutions were obtained ■ 
with relative ease. Once a converged solution is obtained, optimal solutions for other 
cases with differenet inclination changes or different peak heating rate constraints can 
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be obtained using this converged solution as initial guesses. However, it is important 
to point out proper scaling of the defects, constraints and variables are essential to get 
converged solutions. For simulations discussed here, converged solutions were 
obtained by using as little as 50 nodes. However, in some cases, converged solutions 
were obtained using 100 nodes. As far as we know, this may be the first time 
converged soultions were obtained for so many independent variables and nonlinear 
constraint equations. This also illustrates how powerful the nonlinear programming 
code and the collocation and Hermite polynomial technique are. 

Finally, it is important to mention again that aeroassisted orbital transfer introduces a 
strong coupling between the vehicle design and the trajectory design as indicated by 
the simulation data. A trajectory that minimizes fuel mass, without attention to heating , 
may require the vehicle to have heavy thermal protection systems. As shown here, an 
optimal design for the total vehicle weight may be obtained as discussed earlier. 
However, if the aeroassisted transfer is to be prefered to all propulsive transfer, it must 
offer a reduction in fuel mass greater than the increase in thermal protection mass. 

As far as minimum fuel is concerned, the reference cases investigated in Reference 17 
provided more fuel savings as expected. But for the over all trade-off studies, the peak 
heating rate, dynamical pressure, maximum g forces, and fuel mass have to be 
considered. May be, it is also important to point out that the problems investigated 
here is to assume that all plane changes are propulsive and outside the atmosphere 
and that the aeroassisted atmospheric flight is planar. This case is most beneficial 
from the thermal protection point of view and must be considered in the over all trade- 
off studies. 
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Fig. 3 Time History of Altitude 


Fig. 4 Time History of Velocity 
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